# ----------------- #
# Leinden models    #
# ----------------- #

## Clear working space
rm(list = ls())

## Imputed data
war.data <- readRDS("data/war_data_replication.rds")



# Create function to set data in a range of 0 to 1
range01 <- function(x){(x-min(x, na.rm = T))/(max(x, na.rm = T)-min(x, na.rm = T))}

## Load in the dyad data
dyad_data <- readRDS("intl_order/network_layers_replication.rds")

## Make sure the sender and reciever is unique
dyad_data <- distinct(
  dyad_data, 
  ccode1, ccode2, year, .keep_all = T 
) %>% transform(
  mid = cowmidongoing 
) %>% transform(
  fatal_mid = ifelse(mid == 1 & fatality > 0, 1, 0)
) 

## Set missing distance data to missing
dyad_data$contig[dyad_data$mindist == 0] <- 1

cor_data <- dplyr::select(
  dyad_data,
  atop_defense, atop_offense, pta, igo, dcaV2
)

cor_data <- cor_data[complete.cases(cor_data),]
cor_data <- cor(cor_data)
cor_data <- data.frame(cor_data)
cor_data$sender <- row.names(cor_data) 
cor_data <- reshape2::melt(
  cor_data, 
  id = "sender"
)

cor_data <- mutate(
  cor_data, 
  sender = as.character(sender),
  sender = ifelse(sender == "atop_defense", "Defensive\nalliance", sender),
  sender = ifelse(sender == "atop_offense", "Offensive\nalliance", sender),
  sender = ifelse(sender == "dcaV2", "DCA", sender),
  sender = ifelse(sender == "pta", "PTA", sender),
  sender = ifelse(sender == "igo", "IGO", sender),
  
  variable = as.character(variable),
  variable = ifelse(variable == "atop_defense", "Defensive\nalliance", variable),
  variable = ifelse(variable == "atop_offense", "Offensive\nalliance", variable),
  variable = ifelse(variable == "dcaV2", "DCA", variable),
  variable = ifelse(variable == "pta", "PTA", variable),
  variable = ifelse(variable == "igo", "IGO", variable),
  fac_lab = "Layer correlation"
) %>% filter(
  variable != sender
)
unweighted_data <- distinct(
  dyad_data, 
  ccode1, ccode2, year, .keep_all = T
) %>% mutate(
  major_power = cowmaj1  + cowmaj2,
  joint_major_power = if_else(cowmaj1 == 1 & cowmaj2 == 1, 1, 0),
  joint_region = if_else(region1 == region2, 1, 0)
) %>% dplyr::select(
  c(
    ccode1, ccode2, year, contig, 
    pta, bla, dcaV1, dcaV2,
    dipl, igo, 
    atop_defense, atop_offense, atop_neutral, atop_consul,
    joint_region,
    major_power, joint_major_power, cowc1, cowc2,
    fatal_mid, mid, region1, region2
  )
) %>% reshape2::melt(
  id = c(
    "ccode1", "ccode2", "year", "contig", 
    "joint_region", "cowc1", "cowc2",
    "major_power", "joint_major_power", "mid", "fatal_mid", "region1", "region2"
  )
) %>% dplyr::rename(
  layer = variable
) %>% mutate(
  mid_id = paste0(cowc1, "-", cowc2),
  value = if_else(is.na(value), 0, value),
  region_id = paste0(region1, "-", region2)
) 


## Manually fix missing contiguity data
unweighted_data$contig[is.na(unweighted_data$contig) & unweighted_data$ccode1 == 626 & unweighted_data$ccode2 %in% c(625, 530, 500, 501, 484, 490)] <- 1
unweighted_data$contig[is.na(unweighted_data$contig) & unweighted_data$ccode2 == 626 & unweighted_data$ccode1 %in% c(625, 530, 500, 501, 484, 490)] <- 1
unweighted_data$contig[is.na(unweighted_data$contig) & unweighted_data$ccode1 == 626 & !(unweighted_data$ccode2 %in% c(625, 530, 500, 501, 484, 490))] <- 0
unweighted_data$contig[is.na(unweighted_data$contig) & unweighted_data$ccode2 == 626 & !(unweighted_data$ccode1 %in% c(625, 530, 500, 501, 484, 490))] <- 0
unweighted_data$contig[unweighted_data$year >= 2011 & unweighted_data$ccode1 == 625 & unweighted_data$ccode2 %in% c(500, 501, 490)] <- 0
unweighted_data$contig[unweighted_data$year >= 2011 & unweighted_data$ccode2 == 625 & unweighted_data$ccode1 %in% c(500, 501, 490)] <- 0
sum(is.na(unweighted_data$contig))

## Find data to index on
max(subset(unweighted_data, contig == 1)$year)
sum(is.na(unweighted_data$major_power))

## Major power variable creation
unweighted_data$major_power <- factor(
  unweighted_data$major_power,
  levels = c("0", "1", "2")
)

## Create a year and year squared variable
input_data <- unweighted_data
input_data$year_sqr <- input_data$year ^ 2

## Normalize IGO data
igo_temp <- filter(
  input_data, 
  grepl("igo", layer)
) 

range01_data <- list()
for (i in 1:length(unique(igo_temp$layer)))
{
  temp_subset <- filter(
    igo_temp,
    layer == unique(igo_temp$layer)[i]
  )
  for (z in 1:length(unique(temp_subset$year)))
  {
    if (any(subset(temp_subset, year == unique(temp_subset$year)[z])$value != 0))
    {
      subset_by_year <- filter(
        temp_subset, 
        year == unique(temp_subset$year)[z]
      )
      
      ## Round data for replication purposes
      ## This stuff varies slightly by operating system
      subset_by_year <- transform(
        subset_by_year,
        value = round(range01(value), 2)
      )
      range01_data[[length(range01_data) + 1]] <- subset_by_year
    }
  }
}

## Bring back in the normalized IGO data
igo_temp <- do.call(rbind, range01_data)
input_data <- filter(
  input_data,
  !grepl("igo", layer)
) %>% bind_rows(
  igo_temp
) %>% mutate(
  layer = as.character(layer)
)

## Original network layers
input_data <- filter(
  input_data,
  !is.na(value),
  layer %in% c(
    "pta", 
    "atop_offense", 
    "atop_defense", 
    "igo",
    "dcaV2"
  )
) 

gc()

## Loop through the layers to weight data
layer_summary <- list()
for (i in 1:length(unique(input_data$layer)))
{
  
  
  temp_data <- filter(
    input_data, 
    layer == unique(input_data$layer)[i]
  ) %>% distinct(
    ccode1, ccode2, year, .keep_all = T
  ) %>% mutate(
    dyad_id = paste0(cowc1, "-", cowc2)
  )
  
  
  years_unique <- data.table::setDT(temp_data)[,.(n = length(unique(value))), by = 'year']
  temp_data <- filter(
    temp_data,
    year >= min(years_unique$year[years_unique$n > 1])
  )
  
  suppressWarnings({suppressMessages({layer_model <- fixest::feols(value ~ contig | year + ccode1 + ccode2, data = temp_data)})})
  temp_data$pred_values <- predict(layer_model)
  layer_summary[[i]] <- data.frame(temp_data)
}

gc()
## Merge back in weighted data
layer_summary_bind <- do.call(bind_rows, layer_summary)
layer_summary_bind[is.na(layer_summary_bind)] <- 0
layer_summary_bind$pred_values <- round(layer_summary_bind$pred_values, 2)
layer_summary_bind$tie <- layer_summary_bind$value - layer_summary_bind$pred_values


## Normalize the ties by number of layers
weighted_data <- group_by(
  layer_summary_bind,
  year
) %>% mutate(
      n = length(unique(layer))
  ) %>% ungroup(
    .
  ) %>% mutate(
    tie = tie/n,
    tie = round(tie, 2) ## Round for replication purposes across operating systems
  )

## Sum the weighted ties
weighted_data <- data.table::setDT(weighted_data)[,.(tie = sum(tie)), by = 'cowc1,cowc2,year']

## Rename weigthed ties
weighted_data <- dplyr::rename(
  weighted_data,
  namea = cowc1, 
  nameb = cowc2
) %>% data.frame(.)
  
## Create network data for analysis
adjacency_matricies <- list()
for (v in min(unweighted_data$year):2014)
{
  ## Normalize data by year
  temp_edgelist <- filter(
    weighted_data,
    year == v
  )
  
  ## Normalize between 0 and 1
  temp_edgelist$tie[temp_edgelist$tie < 0] <- NA
  temp_edgelist$tie  <- range01(temp_edgelist$tie)
  temp_edgelist$tie <- round(temp_edgelist$tie, 2)
  temp_edgelist <- arrange(
    temp_edgelist,
    namea, nameb
  ) %>% distinct(
    namea, nameb, tie
  )
  
  ## Set negative edge weights to zero
  temp_edgelist$tie[is.na(temp_edgelist$tie)] <- 0
  temp_net <- graph_from_edgelist(as.matrix(temp_edgelist[,c(1:2)]), directed = F)
  E(temp_net)$weight <- temp_edgelist$tie
  
  ## Create an adjacency matrix
  temp_adj <- as_adjacency_matrix(temp_net, attr="weight", sparse = T)
  temp_adj <- as.matrix(temp_adj)
  temp_adj <- temp_adj[order(rownames(temp_adj)),order(colnames(temp_adj))]
  adjacency_matricies[[length(adjacency_matricies) + 1]] <- as.matrix(temp_adj)
}


gc()
set.seed(9732)
## Interlayer community detection
community.blondel <- lapply(adjacency_matricies, function(x){ 
  g <- graph_from_adjacency_matrix(x, mode = "undirected", weighted = T)
  comm <-  cluster_leiden(g, weights = E(g)$weight, n_iterations = 1000, objective_function = "modularity")
  glag.info <- data.frame(
    names = comm$names,
    group = comm$membership
  ) 
    return(glag.info)
}) 
  
## Define year in data and new group names by year
for (i in 1:length(community.blondel))
{
  community.blondel[[i]]$year <- i + min(unweighted_data$year) - 1
  if (i != 1)
  {
    community.blondel[[i]]$group <- community.blondel[[i]]$group + max(community.blondel[[i - 1]]$group) 
  }
}

## Named variables  
community.blondel <- do.call(bind_rows, community.blondel)
community.blondel$group <- as.numeric(as.factor(community.blondel$group))


## Check for missing
missing_data <- unique(
  c(
    paste0(subset(dyad_data, year <= 2014)$cowc1, "-", subset(dyad_data, year <= 2014)$year),
    paste0(subset(dyad_data, year <= 2014)$cowc2, "-", subset(dyad_data, year <= 2014)$year)
  )
)

missing_data[!(missing_data %in% paste0(community.blondel$names, "-", community.blondel$year))]

## Create interlayer network ties
group_matrix <- matrix(0, nrow = length(unique(community.blondel$group)), ncol = length(unique(community.blondel$group)))
colnames(group_matrix) <- rownames(group_matrix) <- unique(community.blondel$group)
    
gc()    
max_year <- plyr::ddply(
    community.blondel, 
    ~group, 
    summarize, 
    max_year = max(year)
  ) %>% arrange(
        max_year
  )

## Create data to fill in 
group_names <- unique(colnames(group_matrix))
edgelist <- data.frame(from = NA, to = NA, weight = NA)
for (i in 1:nrow(max_year))
{
  names_temp <- subset(community.blondel, group == max_year$group[i])$names
  year_temp  <- max_year$max_year[i] + 1
  temp_data <- subset(community.blondel, names %in% names_temp & year == year_temp)
  from <- max_year$group[i]
  to <- unique(temp_data$group)
  weights <- c()
      
  if (length(to) > 0)
  {
    ## Jaccard similarity for each layer
    for (z in 1:length(to))
    {
      union_data <- union(names_temp, subset(temp_data, group == to[z])$names)
      sender_vec <- ifelse(union_data %in% names_temp, 1, 0)
      reciever_vec <- ifelse(union_data %in% subset(temp_data, group == to[z])$names, 1, 0)
      weights[z] <- length(intersect(names_temp, subset(temp_data, group == to[z])$names))/length(union_data)
    }
    
    ## Create Jaccard similarity edge list
    edgelist <- bind_rows(
      edgelist,
      data.frame(
        from = rep(from, length(to)),
        to = to,
        weight = weights
      )
      ) %>% filter(
        !is.na(from)
      )
    }
}

## Round for replication purposes
edgelist$weight <- round(edgelist$weight, 2)

## Fill in an empty matrix for edge weights
for (i in 1:nrow(edgelist))
{
  group_matrix[rownames(group_matrix) == edgelist$from[i],
               colnames(group_matrix) == edgelist$to[i]] <- edgelist$weight[i] 
  
  group_matrix[rownames(group_matrix) == edgelist$to[i], 
                   colnames(group_matrix) == edgelist$from[i]] <- edgelist$weight[i] 
}
    
## Create self loop for 2 times Jaccard similarity
diag(group_matrix) <- 2
group_matrix[group_matrix < 0] <- 0
 
## Run interlayer community detection   
gc()    
g <- graph_from_adjacency_matrix(group_matrix, mode = "undirected", weighted = T)
set.seed(477)
comm <-  cluster_leiden(g, weights = E(g)$weight, n_iterations = 1000, objective_function = "modularity")
between.group <- data.frame(
  group = comm$names,
  interlayer = comm$membership
)
    
## Merge in data
community.blondel <- mutate(
  community.blondel,
    group = as.character(group)
  ) %>% left_join(
    between.group,
      by = c(
        "group"
      )
) 
  
## Find group size for data  
group.size <- mutate(
  community.blondel,
    n = 1
  ) %>% aggregate(
    n ~ year + interlayer, data = ., FUN = sum
)

## Merge in group size    
community.blondel <- left_join(
  community.blondel,
  group.size,
    by = c(
      "year", "interlayer"
    )
)

## Are any missing? 
sum(is.na(community.blondel$interlayer))
missing_data[!(missing_data %in% paste0(community.blondel$names, "-", community.blondel$year))]

## Create merging data
ccode_merge <- data.frame(
  ccode = c(dyad_data$ccode1, dyad_data$ccode2),
  names = c(dyad_data$cowc1, dyad_data$cowc2)
) %>% distinct(
  names, ccode
)

community.blondel <- left_join(
  community.blondel, 
  ccode_merge,
  by = "names"
)

community.blondel.leiden <- dplyr::select(
  community.blondel, 
  c(
    year, ccode, interlayer
  )
)

war.data <- left_join(
  war.data,
  community.blondel.leiden, 
  by = c(
    "ccode1" = "ccode",
    "year"
  ) 
) %>% left_join(
  community.blondel.leiden, 
  by = c(
    "ccode2" = "ccode",
    "year"
  )
) %>% mutate(
  out_group = ifelse(interlayer.x != interlayer.y, 1, 0)
)

bivariate_regression <- fixest::feglm(fatal_mid ~ out_group | ccode1 + ccode2 + dyad + year, data = war.data, family = "binomial")
bivariate_lpm <- fixest::feols(fatal_mid ~ out_group | ccode1 + ccode2  + dyad + year, data = war.data)
full_regression <- fixest::feglm(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio | ccode1 + ccode2  + dyad + year, data = war.data, family = "binomial")
full_lpm <- fixest::feols(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio | ccode1 + ccode2  + dyad + year, data = war.data)

full_regression_mid <- fixest::feglm(mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio | ccode1 + ccode2  + dyad + year, data = war.data, family = "binomial")
full_lpm_mid        <- fixest::feols(mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio | ccode1 + ccode2  + dyad + year, data = war.data)

full_glm        <- glm(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio + contig + poly(year, 2), data = war.data, family = "binomial")
full_lm         <- lm(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio + contig + poly(year, 2), data = war.data)

full_glm_mid    <- glm(mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio + contig + poly(year, 2), data = war.data, family = "binomial")
full_lm_mid     <- lm(mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio + contig + poly(year, 2), data = war.data)

summary(bivariate_regression, cluster = c("year", "ccode1", "ccode2", "nameb"))
summary(bivariate_lpm, cluster = c("year", "dyad", "ccode1", "ccode2"))
summary(full_regression, cluster = c("year", "ccode1", "ccode2", "nameb"))
summary(full_lpm, cluster = c("year", "dyad", "ccode1", "ccode2"))
summary(full_regression_mid, cluster = c("year", "ccode1", "ccode2", "nameb"))
summary(full_lpm_mid, cluster = c("year", "dyad", "ccode1", "ccode2"))
summary(full_glm)
summary(full_lm)
summary(full_glm_mid)
summary(full_lm_mid)




save(full_regression, file = "model_output/leiden_full_regression.rda")
save(bivariate_regression, file = "model_output/leiden_simple_regression.rda")
save(bivariate_lpm, file = "model_output/leiden_simple_lpm_regression.rda")
save(full_lpm, file = "model_output/leiden_full_lpm.rda")
save(full_regression_mid, file = "model_output/leiden_full_regression_mid.rda")
save(full_lpm_mid, file = "model_output/leiden_simple_regression_mid.rda")
save(full_glm, file = "model_output/leiden_full_glm.rda")
save(full_lm, file = "model_output/leiden_full_lm.rda")
save(full_glm_mid, file = "model_output/leiden_full_glm_mid.rda")
save(full_lm_mid, file = "model_output/leiden_full_lm_mid.rda")



## Latent factor model with random year, dyad and node
war.data$out_group          <- as.factor(as.character(war.data$out_group))
war.data$major_power_joint  <- as.factor(as.character(war.data$major_power_joint))
war.data$contig             <- as.factor(as.character(war.data$contig))
war.data$democ_joint        <- as.factor(as.character(war.data$democ_joint))
war.data$autoc_joint        <- as.factor(as.character(war.data$autoc_joint))
war.data$year               <- as.factor(as.character(war.data$year))
war.data$namea              <- as.factor(as.character(war.data$ccode1))
war.data$nameb              <- as.factor(as.character(war.data$ccode2))
war.data$dyad               <- as.factor(as.character(war.data$dyad))
war.data$state_count        <- as.factor(as.character(war.data$state_count))
war.data$year_num <- as.numeric(as.character(war.data$year))
war.data$year_num <- war.data$year_num - mean(war.data$year_num)

set.seed(123)
system.time(multi_factor_model <- brm(data = war.data, 
                                      family = bernoulli,
                                      fatal_mid ~ 1 +
                                        (1|year) + (1|state_count) + 
                                        (1|dyad) + (1|namea) + (1|nameb) + 
                                        out_group + major_power_joint + contig + democ_joint + 
                                        autoc_joint + cinc_ratio + poly(year_num, 2), 
                                      prior = c(set_prior("normal(0, 1)", class = "b"),
                                                set_prior("cauchy(0, 1)", class = "sd")),
                                      chains = 8, iter = 2000, warmup = 1500, cores = 8,
                                      control = list(adapt_delta = 0.95)))
summary(multi_factor_model)
save(multi_factor_model, file = "model_output/leiden_multi_factor_model.RData")
multi_sum <- summary(multi_factor_model)
save(multi_sum, file = "model_output/leiden_multi_factor_summary.RData")

extract.brms <- function(model,
                         ...) {
  s <- model
  coefficient.names <- rownames(s$fixed)
  mean <- s$fixed[,1]
  se   <- s$fixed[,2]
  pval <- ifelse(s$fixed[,3] > 0 | s$fixed[,4] < 0, 0.049, 0.51)
  
  coefficient.names <- gsub("1", "", coefficient.names)
  coefficient.names <- gsub("\\[\\[i\\]\\]", "", coefficient.names)
  keep_vars <- which(grepl("memory|poly", coefficient.names) == F)
  coefficient.names[coefficient.names == "Intercept"] <- "Constant"
  coefficient.names[coefficient.names == "out_group"] <- "Out-cooperative community relation"
  coefficient.names[coefficient.names == "democ_joint"] <- "Joint democracy"
  coefficient.names[coefficient.names == "autoc_joint"] <- "Joint autocracy"
  coefficient.names[coefficient.names == "major_power_joint"] <- "Joint major powers"
  coefficient.names[coefficient.names == "cinc_ratio"] <- "Cinc ratio"
  coefficient.names[coefficient.names == "contig"] <- "Contiguity"
  coefficient.names[coefficient.names == "unaffiliated_joint"] <- "Joint unaffiliated"
  coefficient.names[coefficient.names == "lagged_fatal_mid"] <- "lagged fatal MID"
  tr <- createTexreg(
    coef.names = coefficient.names[keep_vars],
    coef = mean[keep_vars],
    se = se[keep_vars],
    pvalues = pval[keep_vars]
  )
  return(tr)
}


extract.btergm <- function(model,
                           ...) {
  s <- summary(model)
  coefficient.names <- rownames(s)
  mean <- s[,1]
  lb <- s[,2]
  ub <- s[,3]
  coefficient.names <- gsub("edgecov\\.", "", coefficient.names)
  coefficient.names <- gsub("\\[\\[i\\]\\]", "", coefficient.names)
  keep_vars <- which(grepl("memory|poly", coefficient.names) == F)
  coefficient.names[coefficient.names == "edges"] <- "Constant"
  coefficient.names[coefficient.names == "out_group"] <- "Out-cooperative community relation"
  coefficient.names[coefficient.names == "democ_joint"] <- "Joint democracy"
  coefficient.names[coefficient.names == "autoc_joint"] <- "Joint autocracy"
  coefficient.names[coefficient.names == "major_power_joint"] <- "Joint major powers"
  coefficient.names[coefficient.names == "cinc_ratio"] <- "Cinc ratio"
  coefficient.names[coefficient.names == "contig"] <- "Contiguity"
  coefficient.names[coefficient.names == "joint_autoc_sensitivity"] <- "Joint autocracy"
  coefficient.names[coefficient.names == "nodecov.state_count"] <- "Count of states"
  #coefficient.names[coefficient.names == "gwdeg.fixed.1"] <- "Geometrically downweighted degree"
  #coefficient.names[coefficient.names == "gwdeg.fixed.1"] <- "Geometrically downweighted degree"
  coefficient.names[coefficient.names == "joint_democ_sensitivity"] <- "Joint democracy"
  tr <- createTexreg(
    coef.names = coefficient.names[keep_vars],
    coef = mean[keep_vars],
    ci.low = lb[keep_vars],
    ci.up = ub[keep_vars]
  )
  return(tr)
}

extract.fixest <- function(model,
                           include.deviance = TRUE,
                           include.nobs = TRUE,
                           ...) {
  require(fixest)
  s <- summary(model, cluster = c("year", "dyad", "ccode1", "ccode2"))
  coefficient.names <- rownames(s$coeftable)
  keep_vars <- which(grepl("memory|time|poly", coefficient.names) == F)
  coefficient.names[coefficient.names == "edges"] <- "Constant"
  coefficient.names[coefficient.names == "out_group"] <- "Out-cooperative community relation"
  coefficient.names[coefficient.names == "democ_joint"] <- "Joint democracy"
  coefficient.names[coefficient.names == "autoc_joint"] <- "Joint autocracy"
  coefficient.names[coefficient.names == "major_power_joint"] <- "Joint major powers"
  coefficient.names[coefficient.names == "cinc_ratio"] <- "Cinc ratio"
  coefficient.names[coefficient.names == "contig"] <- "Contiguity"
  
  co <- s$coeftable[, 1]
  se <- s$coeftable[, 2]
  pval <- s$coeftable[, 4]
  
  
  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()
  
  if (include.deviance == TRUE) {
    gof <- c(gof, deviance(model))
    gof.names <- c(gof.names, "Deviance")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.nobs == TRUE) {
    n <- s$nobs
    gof <- c(gof, n)
    gof.names <- c(gof.names, "Num. obs.")
    gof.decimal <- c(gof.decimal, FALSE)
  }
  
  
  
  tr <- createTexreg(
    coef.names = coefficient.names[keep_vars],
    coef = co[keep_vars],
    se = se[keep_vars],
    pvalues = pval[keep_vars],
    gof.names = gof.names,
    gof = gof,
    gof.decimal = gof.decimal
  )
  return(tr)
}

extract.glm <- function(model,
                        include.deviance = TRUE,
                        include.nobs = TRUE,
                        ...) {
  s <- summary(model)
  coefficient.names <- rownames(s$coefficients)
  coefficient.names <- rownames(s$coefficients)
  keep_vars <- which(grepl("memory|time|poly", coefficient.names) == F)
  coefficient.names[grepl("intercept", coefficient.names, ignore.case = T)] <- "Constant"
  coefficient.names[coefficient.names == "out_group"] <- "Out-cooperative community relation"
  coefficient.names[coefficient.names == "democ_joint"] <- "Joint democracy"
  coefficient.names[coefficient.names == "autoc_joint"] <- "Joint autocracy"
  coefficient.names[coefficient.names == "major_power_joint"] <- "Joint major powers"
  coefficient.names[coefficient.names == "cinc_ratio"] <- "Cinc ratio"
  coefficient.names[coefficient.names == "contig"] <- "Contiguity"
  
  co <- s$coefficients[, 1]
  se <- s$coefficients[, 2]
  pval <- s$coefficients[, 4]
  
  
  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()
  
  if (include.deviance == TRUE) {
    gof <- c(gof, deviance(model))
    gof.names <- c(gof.names, "Deviance")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.nobs == TRUE) {
    n <- s$df[2]
    gof <- c(gof, n)
    gof.names <- c(gof.names, "Num. obs.")
    gof.decimal <- c(gof.decimal, FALSE)
  }
  
  
  
  tr <- createTexreg(
    coef.names = coefficient.names[keep_vars],
    coef = co[keep_vars],
    se = se[keep_vars],
    pvalues = pval[keep_vars],
    gof.names = gof.names,
    gof = gof,
    gof.decimal = gof.decimal
  )
  return(tr)
}

## hierarchical coefficients
extract.re.model <- function(model,
                             include.loglik = TRUE,
                             include.deviance = TRUE,
                             include.nobs = TRUE,
                             ...) {
  s <- summary(model)
  coefficient.names <- rownames(s$coefficients)
  coefficient.names <- rownames(s$coefficients)
  coefficient.names[coefficient.names == "prop_autoc_mimic"] <- "Autocratic exposure"
  coefficient.names[coefficient.names == "prop_democ_mimic"] <- "Democratic exposure"
  coefficient.names[coefficient.names == "prop_exec_constraint"] <- "Executive constraint exposure"
  coefficient.names[coefficient.names == "prop_religious_freedom_mimic"] <- "Religious freedom exposure"
  coefficient.names[coefficient.names == "prop_rest_pol_part"] <- "Restricted political participation exposure"
  coefficient.names[coefficient.names == "prop_gov_censorship_mimic"] <- "Free media exposure"
  
  
  
  coefficient.names[coefficient.names == "autocratic_neighbor"] <- "Autocratic neighbor (%)"
  coefficient.names[coefficient.names == "exec_constraint_neighbor"] <- "Executive constraint neighbor  (%)"
  coefficient.names[coefficient.names == "democratic_neighbor"] <- "Democratic neighbor  (%)"
  coefficient.names[coefficient.names == "religious_freedom_neighbor"] <- "Religious freedom neighbor (%)"
  coefficient.names[coefficient.names == "rest_pol_part_neighbor"] <- "Restricted political participation neighbor (%)"
  coefficient.names[coefficient.names == "gov_censorship_neighbor"] <- "Free media neighbor (%)"
  
  coefficient.names[coefficient.names == "(Intercept)"] <- "Constant"
  
  coefficient.names[coefficient.names == "civil_war"] <- "Civil war"
  coefficient.names[coefficient.names == "log(gle_cgdpc)"] <- "log(GDPC)"
  coefficient.names[coefficient.names == "wdi_gdpcapgr"] <- "Annual growth GDPC"
  coefficient.names[coefficient.names == "unaffiliated"] <- "Isolate"
  
  co <- s$coefficients[, 1]
  se <- s$coefficients[, 2]
  pval <- s$coefficients[, 4]
  
  
  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()
  if (include.loglik == TRUE) {
    lik <- logLik(model)
    gof <- c(gof, lik)
    gof.names <- c(gof.names, "Log Likelihood")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.deviance == TRUE) {
    gof <- c(gof, deviance(model))
    gof.names <- c(gof.names, "Deviance")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.nobs == TRUE) {
    n <- nrow(model@frame)
    gof <- c(gof, n)
    gof.names <- c(gof.names, "Num. obs.")
    gof.decimal <- c(gof.decimal, FALSE)
  }
  
  
  
  tr <- createTexreg( 
    coef.names = coefficient.names[!grepl("poly", coefficient.names)], 
    coef = co[!grepl("poly", coefficient.names)],
    se = se[!grepl("poly", coefficient.names)],
    pvalues = pval[!grepl("poly", coefficient.names)],
    gof.names = gof.names,
    gof = gof,
    gof.decimal = gof.decimal
  )
  
  return(tr)
}


load("model_output/leiden_full_regression.rda")
load("model_output/leiden_simple_regression.rda")
load("model_output/leiden_simple_lpm_regression.rda")
load("model_output/leiden_full_lpm.rda")
load("model_output/leiden_full_regression_mid.rda")
load("model_output/leiden_simple_regression_mid.rda")
load("model_output/leiden_full_glm.rda")
load("model_output/leiden_full_lm.rda")
load("model_output/leiden_full_glm_mid.rda")
load("model_output/leiden_full_lm_mid.rda")
load("model_output/leiden_multi_factor_summary.RData")

library(texreg)
full_fe_logit <- extract.fixest(full_regression)
full_fe_lpm   <- extract.fixest(full_lpm)
full_glm      <- extract.glm(full_glm)
full_lm       <- extract.glm(full_lm)
full_lfm      <- extract.brms(multi_sum)

mid_full_fe_logit <- extract.fixest(full_regression_mid)
mid_full_fe_lpm   <- extract.fixest(full_lpm_mid)
mid_full_glm      <- extract.glm(full_glm_mid)
mid_full_lm       <- extract.glm(full_lm_mid)


si_table_9 <- texreg(list(full_lm, full_glm, full_fe_lpm, full_fe_logit, full_lfm, 
            mid_full_lm, mid_full_glm, mid_full_fe_lpm, mid_full_fe_logit), stars = 0.05, digits = 2,
            reorder.coef = c(2:7, 1))


write.table(si_table_9, file = "data_visualization/si_table_9.txt")
